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We study leasing by voids in Cubic Galileon and Nonlocal gravity cosmologies, which are examples of 
theories of gravity that modify the leasing potential. We find voids in the dark matter and halo density fields of 
N-body simulations and compute their leasing signal analytically from the void density profiles, which we show 
are well fit by a simple analytical formula. In the Cubic Galileon model, the modifications to gravity inside voids 
are not screened and they approximately double the size of the leasing effects compared to GR. The difference 
is largely determined by the direct effects of the fifth force on leasing and less so by the modified density 
profiles. For this model, we also discuss the subtle impact on the force and leasing calculations caused by the 
screening effects of haloes that exist in and around voids. In the Nonlocal model, the impact of the modified 
density profiles and the direct modifications to leasing are comparable, but they boost the leasing signal by only 
~ 10%, compared with that of GR. Overall, our results suggest that leasing by voids is a promising tool to test 
models of gravity that modify leasing. 


I. INTRODUCTION 

Despite the success of General Relativy (GR) in passing all 
currently available solar system tests of gravity m, there is 
growing interest in the theoretical lIHO and observational [I4F 
121 aspects of theories beyond GR. There are two main reasons 
for this. Firstly, the simple fact that GR has not been tested 
on scales larger than the solar system means that, in fact, one 
makes a huge extrapolation of the regime of validity of the 
theory when one uses it (as it is common) in cosmological 
studies. The gravitational law should, therefore, be put to test 
on larger scales, and modified gravity models help to identify 
the types of imprints that modifications to gravity can leave 
on observables. Secondly, there is currently no theoretically 
appealing explanation for the nature of the dark energy that 
is responsible for the accelerated expansion of the Universe. 
In the standard A-Cold Dark Matter (ACDM) cosmological 
model, the role of the dark energy is attributed to a simple cos¬ 
mological constant A, but the smallness of its value remains a 
mystery. Models of modified gravity can explain the acceler¬ 
ation without A, thereby providing extra motivation for their 
study. 

The majority of modified gravity models predict the exis¬ 
tence of extra degrees of freedom (often of the scalar type) that 
fifth forces felt by the matter fields. Consequently, a 
major difficulty in building models of modified gravity comes 
from making them compatible with the stringent solar system 
bounds. The latter constrain the fifth force to be extremelly 
small, and hence, cosmologically uninteresting. A popular 
way out of this relies in building models that possess what are 
commonly referred to as screening mechanisms. In short, the 
idea is to construct models where the equations of the scalar 
field become highly nonlinear in regions of high density (like 
the solar system). The presence of the nonlinearities acts to 
suppress the magnitude of the fifth force. On larger scales, 
where the density is low, the fifth force effects become man¬ 
ifest and potentially detectable. On these large scales, the 
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scalar field equation can be linearized to look like a Pois¬ 
son equation. Examples of screening mechanisms include the 
chameleon mechanism ||2l which operates in the popular Hu- 
Sawicki f{R) lEl gravity model; the Vainshtein mechanism 
iHni which operates in the Dvali-Gabadadze-Porrati (DGP) 
m and Galileon ifTMBI models; the K-moufiage screening 
II1614201 and disformal screening 112111221 . 

Due to the suppression effects of the screening, it is best 
to devise observational tests that focus on large scales or low- 
density regions, where the screening is less efficient 1231 l24l . 
For instance, recent studies have shown that the amplitude 
of the cosmic microwave background (CMB) lensing poten¬ 
tial l20l I25lj77ll and cosmic shear 1^ power spectra 
(which are sensitive to the projected matter distribution on 
large scales) are, indeed, a sensitive probe of modified grav¬ 
ity. The cross-correlation of galaxy positions with the lensing 
shear field can also help to constrain modified gravity oniiiii. 
The integrated Sachs-Wolfe (ISW) effect, which probes the 
time variation of large scale gravitational potentials, consti¬ 
tutes another good example of constraining gravity away from 
the regimes where the screening is at play l25ll26ll32H36l . On 
the other hand, although the amplitude of the matter power 
spectrum on large scales is also affected by the modifications 
to gravity, the uncertainties about galaxy bias undermine the 
possibility of obtaining tight constraints (see e.g. Sec. IV. D 
of Ref. EtI for a discussion). On mildly nonlinear scales 
(2 — 20Mpc), several recent studies have found that the pe¬ 
culiar velocities of galaxies are also very sensitive to the pres¬ 
ence of fifth forces lISSlBTl . These scales are typically associ¬ 
ated with the infall regions of massive galaxy clusters, which 
are located sufficiently far away from the cluster center for 
the screening to have a smaller impact. In general, inside the 
virial radius of galaxy clusters (< IMpc), it becomes harder 
to find the effects of the fifth force (see e.g. Ref. Il42l l. 

Here, we focus on cosmic voids, which are the regions of 
the Universe where the density is the lowest, and hence, where 
one expects fifth force effects to be maximal. Despite being 
potentially good probes of gravity, voids have only recently 
become the object of dedicated studies in modified gravity 
Il43l - l4^ . In particular. Ref. Il47l showed that the lensing signal 
from voids in f{R) gravity is modified relative to ACDM, via 
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the modifications induced by the fifth force to the void density 
profiles. This result is particularly timely as it can be linked 
to the recent work of Refs. EHEQl, who have independently 
detected the lensing signal associated with cosmic voids. This 
therefore opens the prospect of developing new tests of grav¬ 
ity using the lensing signal in and around voids. In terms of 
lensing, f{R) models (and scalar-tensor theories in general) 
are special in the sense that they do not modify the lensing 
signal directly. In these models, the amplitude of the fifth 
force vanishes for relativistic particles like photons. In other 
words, any modifications to lensing arise through changes in 
the mass distribution, and not due to changes to the photon 
geodesic equation. For this reason, one expects that lensing by 
voids can serve as a stronger probe of models that also mod¬ 
ify the photon geodesic equation. Examples of such models 
include Nonlocal gravity EBISSI, Galileon gravity GMia, 
massive gravity EMI, K-mouflage gravity ifTMT^IIOll. Ki¬ 
netic Gravity Braiding llTSl l64l l65l and several other special 
cases of Horndeski’s general model 1661 . One of our goals 
here is to investigate the lensing signal from voids in some of 
these models. 

As working cases, we focus on the Cubic Galileon model 
Ea and the Nonlocal gravity model of Ref. ll5Tl . We make 
use of the N-body simulations performed for these two mod¬ 
els in Refs. 15^l68l . We find voids in the simulations using 
a watershed based algorithm 1^ and investigate the effects 
of the fifth force on the number of voids and on their den¬ 
sity and force profiles. We also put forward a simple fitting 
formula that matches very well the void profiles found in the 
simulations for different variants of the modified gravity mod¬ 
els, for different density tracer types (dark matter and haloes) 
and for a wide range of void sizes. Our formula is a simple 
extension of others used previously GQlEIl, and by having 
more parameters it provides a better fit to our simulation re¬ 
sults. The formula admits a closed expression (in terms of 
hypergeometric functions) for the mass within a given radius, 
which makes it convenient to use in force profile calculations 
and lensing studies. When we assess the impact of the fifth 
force on the lensing signal, we take into account its effect on 
both the void density profiles and the calculation of the lens¬ 
ing observables themselves. Our goal is to provide intuition 
about the potential of lensing by voids to test gravity outside 
the solar system. We do not attempt to make any observation- 
ally conclusive statement, but we do comment on a number of 
extra steps that need to be taken to compare our results with 
observations. 

This paper is organized as follows. In Sec. we intro¬ 
duce the force equations in Nonlocal and Galileon gravity, 
discussing some of their phenomenology. In Sec. [I^ we de¬ 
scribe our N-body simulations and the void finding algorithm, 
and study the effects of the fifth force on the abundance, den¬ 
sity profiles and force profiles of the voids. In Sec. we 
describe the calculation of the lensing signal, and then com¬ 
pute it for the voids found in the simulations. In the same 
section, we also link our findings to recent observational re¬ 
sults, and provide a quick guideline of the steps needed for 
more elaborate comparisons to observations. We summarize 
our results in Sec.lVl 


II. THE MODELS OE GRAVITY 

In this section, we briefly introduce the models of gravity 
that we consider and present the relevant force equations that 
are needed to compute their lensing signal. In the equations 
below we always assume spherical symmetry and work with a 
perturbed Friedmann-Robertson-Walker (FRW) spacetime in 
the Newtonian gauge 

ds^ ={l + 2T'/c2) c^df - (1 - 2$/c 2) dx^, (1) 

where a = 1/(1 -I- z) is the cosmological scale factor (z is the 
redshift) and c is the speed of light. 


A. Nonlocal gravity 

We consider the Nonlocal gravity model of Refs. EllHI- 
Its action is given by 


S = 


IQttG 


J 


- ^TZa-^Tl - 
6 


( 2 ) 


which can be cast in a more familiar (local) form given by 

EMI 


A= 


1 

IGttG 



771^ 

n - ns - a (au + n) 

6 

6 (05 + u)- £^], 


(3) 


where n is the Ricci scalar, g is the determinant of the met¬ 
ric G is Newton’s gravitational constant. Cm is the mat¬ 
ter Fagrangian density, and ^2 are Lagrange multipliers, 
U = —'G\~^n and S = 'G\~'^n are two auxiliary scalar fields 
and □ = is the d’Alembert operator, with Greek in¬ 

dices running over 0, 1, 2, 3. Here, we do not present a de¬ 
tailed discussion about the theoretical aspects of the above two 
actions, but simply caution that their solutions are not com¬ 
pletely equivalent and that care must be taken before matching 
them (see, e.g. Refs. ElElEl for a discussion). 

On the scales relevant for large scale structure formation 
and in the absence of anisotropic stress, the two Newtonian 
potentials are the same (dr = $) and the modifed Poisson 
equation can be written as El|56l 

^ m = 47rGefiPm<5(.R), (4) 

where pm is the cosmological background value of the physi¬ 
cal matter density pm, S = Pmipm — 1 is the density contrast 
and denotes partial differentiation w.r.t. the radial coordi¬ 
nate R. The above equation has the same form as in GR but 
with an effective time-dependent gravitational strength given 
by 


Geff = G 


1 - 


^S{z) 


-1 


> 1 , 


(5) 
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where S is the background part of the field S. The time evolu¬ 
tion of the background quantities in the Nonlocal model have 
to be obtained numerically by integrating the background dif¬ 
ferential equations (see e.g. Refs Il54ll5^ '). The parameter m 
in Eqs. Q and ([^ is controlled by the amount of dark energy 
in the Universe, i.e., in a flat Universe, its value is determined 
by the energy densities of the remaining matter species. This 
means that this Nonlocal gravity model has the same number 
of free parameters as ACDM. For reference, for the model 
parameters used in Ref. Il5^ (e.g. ilmo = 0.30), one has that 
Ges{z = 0)/G« 1.06. 

The fact that this model is characterized by an enhanced 
gravitational strength on all length scales leads to the ques¬ 
tion of whether or not this model is capable of passing solar 
system constraints IT]. In Ref ll5^ . we showed that if the 
gravitational strength of Eq. 0 is used in solar system tests, 
then the model predicts values for the rate of change of the 
gravitational strength, Geff, that are incompatible with cur¬ 
rent lunar laser ranging experiments lISOl . However, the time 
evolution of Geff/G is controled by the background part of 
the held S. This means that if one describes the spacetime 
around the Sun as perturbed Minkowskii (instead of FRW), 
then 5 = 0 —> Geff = 0 , rendering the model compatible 
with current bounds ESI ED. Here, we focus on void size 
scales, which are sensitive to the background expansion, and 
as a result, we use the gravitational strength of Eq. Q when 
computing the model predictions. 


B. Cubic Galileon gravity 

We focus on the Cubic sector of the covariant Galileon 
gravity model ifTTHTSl l82] whose action is given by 

1 ^ - 2C2C2 - -C3C3 - Cm , 

( 6 ) 

where C 2 and C 3 are dimensionless constants, and £2 and £3 
are given by 

£2 = £3 = (7) 

in which (p is the Galileon held, = Mp\Hq, Mpj = 
l/( 87 rG) is the reduced Planck mass squared and Hq — 
lOOh km/s/Mpc is the present-day Hubble expansion rate. In 
flat spacetime, the above action is invariant under the Galilean 
shift dp —)■ dp + (where is a constant four-vector). Fol¬ 
lowing the derivation of Refs. Il68l ETl . the force law in the 
Cubic Galileon model is given by 

G5M{< R) C 3 

TT = —- —■ ® 

where Sp is the spatial perturbation of the Galileon held, p{z) 
is its backround part, and 5M{< R) = -iirpm /q^ 5{x)x^(1x is 
the mass perturbation enclosed in a sphere of radius R. Com¬ 
pared to GR, Eq. 0 has an extra term, which is proportional 
to 







FIG. 1. Representative force profiles in the Cubic Galileon model. 
The top panel shows three example void density profiles. The 
three lowest panels, from top to bottom, show the radial profiles 
of {r./Rf + 1, and for the density profiles shown in 

the top panel. The colors indicate which prediction is associated 
with which density profile. In the bottom two panels, the dashed 
and solid lines correspond, respectively, to the predictions of the lin¬ 
earized (cf. Eq. regime (ii)) and full (cf. Eqs. {nj and GD) 
Cubic Galileon models, as labelled. The meaning of a negative am¬ 
plitude for the force is that it points outwards. The units in the bottom 
two panels are (km/s)^/i^Mpc“^. 
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with 


16 ^pi 

T I3iI32M^ 


G5M{< R), 


where /3i and (32 are two dimensionless functions of time. The 
quantity r* is a radial scale, which is often referred to as the 
(10) Vainshtein radius. From Eqs. ([^, (|^, one can write 
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4 C3 tjp' / 

3MpiM3 hVj 



G6M{< R) 
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( 11 ) 
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( 12 ) 


Contrary to the case of Nonlocal gravity, the Galileon 
model admits analytical solutions for the time evolution of the 
background quantities 1261. The time evolution of the Hubble 
parameter, (p, jSi and /I 2 are given, respectively, by 


ii If r* > 0 and r,/i? ^ 1, then 



/3i 


/32 


Hi 
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iHl/H. 
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6 c 3 

2M^Mpi 


^ ^ 4 " 4(1 — ^mo) 

{(p + 2H0) 


C 2 


4c3 
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/3i. 


(13) 

(14) 
,(15) 

(16) 


As in Ref. Il26l . we take C 2 = — 1 and the other two Galileon 
parameters are determined by flmo as 


e = \/ 6 (l-H^o), (17) 

C3 = 1/(60- (18) 

We take Hmo = 0.28, which is the value used in the simula¬ 
tions of Ref. Il 68 ll . 

From Eq. o, it is possible to identify three regimes for 
the amplitude of the total force in the Cubic Galileon model 
that are relevant for our analysis; 

i In the regime where r* > 0 and r*/i? ^ 1, one has 


and Eq. GD becomes 




Geff = G 1 - - 


C 3 




3 MpiM3 /32 


> 1 


( 22 ) 

(23) 


(/32 < 0 II 68 II I. In this linear regime, which occurs suf- 
hciently far away from massive objects, the force law 
is as in GR, but with an enhanced time-dependent grav¬ 
itational strength. This is similar to the force law of 
the Nonlocal gravity model, albeit with a different time 
evolution for Geff. In particular, in the Galileon model, 
and for the model parameters used in Ref. Il68l . one has 
Geff (2 = 0)/G ~ 2, which is subtantially stronger than 
the « 6 % enhancement in the Nonlocal model. 


iii Einally, there is a third regime characterized by r* < 0 
and |r*/i?| ~ G{\). In this regime, which occurs when¬ 
ever the mass perturbation becomes negative (as it does in 
voids), the total force can be written as 



3 




<1, (19) 


and, as a result, Eq. GD can be approximated as 


^ _ GSM{< R) 

R ^ R^ ' ^ ’ 

That is, close to very massive objects (small R and/or 
large mass perturbations, r* cx SM), the force law in the 
Galileon model becomes the same as in GR. This illus¬ 
trates the implementation of the Vainshtein screening ef¬ 
fect that allows this model to satisfy solar system tests of 
gravity. 


^jR 

R 


Gvoid7?) 


SMi< R) 

^3 ’ 


(24) 


where Gvoid(^, 7?) is a time and scale dependent effec¬ 
tive gravitational strength (simply the term between {} in 
Eq. GD), which is larger in magnitude than the gravita¬ 
tional strength of regime (ii), i.e., Gvoid > Ggff. This can 
be checked by noting that 


Q) 


^ + 1-1 


> 1/2 


(25) 


in Eq. ( [TT] l, when < 0 (cf. Eq. (|2T])). Note that in our 
notation, when SM < 0, then the force becomes negative. 
This means that the force points outwards. 
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To help understand the behavior of the fifth force in the 
Galileon model, we show in Fig. [T] the radial profiles of 
{r^/Rf + 1, and /ii, for each 

of the density profiles depicted in the top panel. The density 
profiles are computed using the formula 


S{R' = R/R,) = 5, ^ (26) 

1 + {R'/S2f 


where Ry is the void radius (whose exact value is not im¬ 
portant for the discussion here) and Sy, a, 13, si and S 2 are 
fitting parameters. Figure[^in the Appendix shows the impact 
that each of the five parameters of Eq. (j2^ has on the density 
profiles (and on the associated lensing signal, whose calcula¬ 
tion is explained in Sec. |IV A| ). In Sec. |IIIC] we shall see that 
this formula provides a very good fit to the void density pro¬ 
files found in the N-body simulations. The mass perturbation, 
dM{< R) = 4Trpm 3(x)x^dx, admits a closed formula 
given by 


SM{< R) = dirpm 


i?3 

3(a-f 3) 


{a + 3) 2 F 1 


3 /3 + 3 
’/3’ /3 


where 2 F 1 is the Gauss hypergeometric series function. This 
formula for 5M{< R) facilitates straightforward calculation 
of the force profiles. In the bottom two panels of Fi g. [T] 
the solid curves show the result obtained by using Eqs. ( |1 l| l 
and ( [T2] i (which we call the full solution), whereas the dashed 
curves show the result associated with the regime (ii) above 
(Eq. ( |22] i, which we call the linear solution). Eor the cases 
shown, for i?' > 1.2 (r* > 0), the full and the linear so¬ 
lutions for the total force, $,/j, are roughly the same, which 
is as expected since these radial scales correspond to the lin¬ 
ear regime (ii) discussed above. A more careful inspection of 
those scales shows that the solid curves underpredict slightly 
the dashed ones. This is due to the Vainshtein screening mech¬ 
anism, which acts to somewhat suppress the full solution. On 
the other hand, for R' < 1.2 (r* < 0), the full solution en¬ 
ters regime (iii), and as expected, the force becomes larger in 
magnitude (more negative) compared to the linear case. 

There is one peculiar aspect about regime (iii) that is worth 
discussing with more detail. The quantity (r»/i?)^ + 1 ap¬ 
pears in Eqs. (Ill and ( [T^ as the argument of square-roots. 
This implies that the amplitude of the force becomes a com¬ 
plex number whenever (r*/i?)^ -f 1 < 0, which is not a phys¬ 
ical result. This problem has been already encountered in the 
N-body simulations of Ref. Il68l . where the authors circunvent 
the absence of real solutions by adopting the ad-hoc fix of set¬ 
ting {r^f/RY -(-1 = 0 whenever it becomes negative. Here, we 
shall implement the same procedure, as we wish to compare 
some of our results to those of Ref. IMl . In Eig.[2 the imple¬ 
mentation of this "fix" is noted by the kinks in and spikes 
in V^<h, for the solid green and solid red curves. The density 
inside the void depicted by the blue line is not low enough for 
{r^/RY -I- 1 to cross zero, and hence, the problem is not seen. 


It is important to try to understand the implications of 
the existence of complex solutions for the fifth force in the 
Galileon model. Here, we note that Eqs. GD and ( [T^ are ob¬ 
tained under the approximations that the perturbed fields are 
weak and quasi-static. In the weak field approximation, one 
neglects terms that involve the perturbed fields and their first 
spatial derivatives, over their second derivatives. The quasi¬ 




f o;-l-3 cy j3 d 



(27) 


static approximation amounts to neglecting the time varia¬ 
tion of the perturbed quantitie^ e.g., (p{a,R) = <^(a) -I- 
6(p{a,R) « <p{a). In the very low-density regions that char¬ 
acterize voids, one expects the weak-field approximation to 
still hold, but the quasi-static one may not (see Refs. Il84ll85ll 
for discussions about the quasi-static limit, and Refs. Il86ll87ll 
for work beyond this in N-body simulations). One may specu¬ 
late that the terms which are neglected in the quasi-static limit 
are actually responsible for keeping the fifth force real for 
all density values when they are present. Interestingly, how¬ 
ever, the recent work of Ref. Il88l has shown that the problem 
remains even after relaxing the quasi-static approximation. 
In particular, the authors find that the time derivative of the 
Galileon field perturbation becomes singular when the quasi¬ 
static solution becomes a complex number (see Ref. Il88l for 
the details). This suggests that the breakdown of the quasi¬ 
static solutions may well be associated with a true instabil¬ 
ity of the Cubic Galileon model. Moreover, the earlier work 
of Ref. Il89l also unveiled some instabilities in general Vain¬ 
shtein solutions around static spherically symmetric sources 
(see Ref. Il89l for the details). Here, we shall keep these dis¬ 
cussions in mind but proceed by retaining the ad hoc fix of 
Ref 1681 . Our treatment of the equations of Cubic Galileon 
can be viewed as a toy model that we use to illustrate the ef¬ 
fects of modified gravity in the properties of voids. 


III. VOIDS IN T HE SIMULATIONS 

A. Outline of the simulations and void finder 

We make use of the N-body simulations of the Cubic 
Galileon and Nonlocal gravity models presented in Refs. l56] 
l68l . Eor the case of the Nonlocal model, the simulations 


' These approximations are also used in the Nonlocal gravity model to obtain 
Eq. |4}, but in this model the problem of complex solutions does not aiise. 
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TABLE I. Summary of the model variants of the Cubic Galileon and 
Nonlocal models studied in this paper. 


Model 

Expansion history 

Force law 

Full Galileon 
Linear Galileon 

Galileon 

Galileon 

Galileon 

Eqs. (IT] 

and 1 12 

22| 

Full Nonlocal 

Nonlocal 

Nonlocal 

Eq. 

G1 



were run with a modified version of the Adaptive Mesh Re- 
hnement (AMR) and publicly-available RAMSES code Il90l . 
The modihcations involved (i) changing the code to interpo¬ 
late the expansion rate in the Nonlocal model using a table; 
(ii) rescaling the gravitational force computed in the default 
code by the value of Geff/G in Eq. (j^, which is also inter¬ 
polated from tabulated values. For the Cubic Galileon model, 
the simulations were performed with the ECOSMOG code ED, 
which is also based on RAMSES. The ECOSMOG code con¬ 
tains additional subroutines that solve the equation of motion 
of the Galileon scalar held via Newton-Gauss-Seidel itera¬ 
tive relaxations on the AMR grid. Once the values of the 
Galileon held are found on the grid, its spatial gradient is 
obtained by hnite differencing to determine the hfth force. 
The modihed background expansion history is also consis¬ 
tently taken into account by the ECOSMOG code. We refer the 
reader to Refs. ED 1211 for more details about how to solve 
Galileon-like equations in adaptively rehned grids (see also 

Refs, limii). 

Hereonin, we shall analyse the results from three variants 
of the Cubic Galileon model and two variants of the Nonlocal 
model, as listed in Table |I] For the Galileon, we call these the 
full, linear and QCDM variants. The QCDM model is char¬ 
acterized by having the background expansion of the Galileon 
model, but the gravitational law of GR. The linear model is the 
same as QCDM, but with the effective gravitational strength 
Geff of Eq. Finally, the full model is, as the name sug¬ 
gests, the Galileon model with its modihed background and 
scale-dependent (with screening) hfth force. Comparing the 
results of the full and linear variants allows one to measure 
the impact of the scale dependence of the hfth force, while 
the QCDM model serves as the reference against which one 
can measure the effects of the modihed gravitational law[^ 
Similarly, for the Nonlocal model, we also have the equiva¬ 
lent QCDM and full model variants. For the Nonlocal model, 
as there is no screening, there is no distinction between the 
linear and the full models. 

We show results from simulation boxes of side 400Mpc//i 
for the Galileon, and 200Mpc/h for the Nonlocal model, both 


^ We do not use standard ACDM as the reference model since the latter dif¬ 
fers from the Galileon model also in the time evolution of the cosmological 
background. Here, we are interested on the effects of the fifth forces alone, 
which is why we use the QCDM variant. 


with 512^ dark matter tracer particles (these were the boxes 
used in Refs. Il5^l68l f. Each of the model variants was simu¬ 
lated hve times using different realizations of the initial den¬ 
sity held. We use the variance across the realizations to com¬ 
pute errorbars. When hnding voids in the simulations, we 
shall also make use of DM haloes found in the simulations. 
Our halo catalogues were obtained with the publicly avail¬ 
able Rockstar code 1951 . which is a phase space friends-of- 
friends based halo hnder. The number density of the haloes 
we consider is nhaio = 5 x 10“^/i^/Mpc^ and rihaio = 
5 X 10“^h^/Mpc^ for the Galileon and Nonlocal gravity sim¬ 
ulations, respectively. This is roughly the number density of 
haloes after retaining only those haloes whose mass is at least 
100 times the particle mass, Mp. This minimum halo mass 
is lOOMp « 4 X W^Mq/Ii and lOOMp « 5 x W'^Mq/Ii, 
for the simulations of the Cubic Galileon and Nonlocal mod¬ 
els, respectively]^ We refer the reader to Refs. l37l 13^IMl 
for further details about the properties of dark matter haloes 
in these models. 

We hnd voids using the Watershed Void Finder (WVF) 
method of Ref. Il69]l . Our code takes as input the dis¬ 
crete tracer distribution, which in our case are DM particles 
and/or DM haloes, to construct a continuous volume-weighted 
density held using a Delaunay Tessellation Field Estimator 
(DTFE) method Il96ll97ll . For computational convenience, the 
DTFE held is sampled onto a regular grid, whose cell size 
is of the order of the mean distance between tracers. The 
grid density held is smoothed with a Gaussian hlter of size 
2 Mpc//i to reduce small scale features that could lead to spu¬ 
rious voids [69]. In the language of the watershed technique, 
the resulting density held is viewed as a landscape that will be 
Hooded by a rising level of water. The regions around every 
local minima of the density held are called catchment basins 
(where water collects) and will be identihed as the voids. As 
the water level rises, the basins grow and, eventually, neigh¬ 
bouring basins meet at the higher-density ridges that separate 
them. These ridges mark the boundary of each basin/void, 
and are associated with the hlaments and walls of the cos¬ 
mic web ll98l l99l . The process stops when the water level 
reaches the global maximum of the density held, by the end 
of which all basin/void boundaries have been identihed. To 
overcome watershed over-segmentatiorj^ ridges whose den¬ 
sity constrast is (5 < —0.8 are not classihed as void bound¬ 
aries, as such low density boundaries are indicative of sub¬ 
voids that have merged EUdTO). An appealing aspect of 
the watershed method is that it makes no a priori assump¬ 
tions on the size, shape or mean underdensity of the voids (see 
Ref. IIQII for a comparison study of different void finders). 

As is customary in void studies, we define the effective void 
radius Ry as the radius of a sphere whose volume is the same 


^ Note that due to the different growth of structure, the halo mass function 
differs between the different variants of the models. The halo catalogues of 
the different variants were cut at slightly different mass values to yield the 
same number density of haloes. 

This refers to avoid finding too many small voids inside a large underdense 
region, where in fact the whole underdense region should be classified as a 
single void that resulted from the merging of smaller ones. 
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as the volume of the watershed void. We take the center of 
the void to be the location of the barycenter which we define 
as Tbarycenter = where f* is the position of each 

grid cell identified as part of the void and A^ceii is the total 
number of grid cells associated with void. We evaluate the 
density profiles of the voids using the DM density field (for 
voids found in both the DM and halo density fields) since this 
is the mass distribution that determines the lensing signal. In 
what follows we limit ourselves to analysing the simulation 
results at z = 0. 


B. Void size function 

Figure shows the cumulative size function of the voids 
found in the simulations of the Cubic Galileon (left panels) 
and Nonlocal (right panels) gravity models. For both models, 
the void population depends on the tracer type used. In partic¬ 
ular, DM density field voids (circles) are smaller and, in total, 
are found in greater number than voids in the halo density 
field (squares). This follows straightforwardly from the fact 
that the distribution of collapsed haloes is sparser than that of 
the DM particles. It is also noteworthy that, for the same type 
of tracer, we find larger voids in the Cubic Galileon than in the 
Nonlocal gravity model. Part of this result is due to the fact 
that the box size used in the simulations of the Galileon model 
(400Mpc//i) is larger than that used in the simulations of the 
Nonlocal model (200Mpc//i). One should therefore bear this 
difference in the box size in mind when comparing the results 
between the two gravity models. 

In terms of the relative difference to QCDM, the full and 
linear variants of the Galileon model predict an enhance¬ 
ment of the order 10% — 20% for the larger DM field voids 
(15Mpc//i ^ Rv ^ 20Mpc//i). This is due to the enhanced 
gravity of these models which boosts the evacuation of matter 
from inside the voids and the formation of large scale struc¬ 
tures. In other words, voids expand faster in the full and linear 
variants, which is why large voids are more abundant. By the 
same reasoning, one should also expect the number of smaller 
voids to be suppressed in the linear and full variants, com¬ 
pared to QCDM. This is because the faster expansion of the 
voids makes it more likely for small neighbouring voids to 
merge into larger ones. In Fig.|^ this suppresion can be seen 
for Ry < 10Mpc//i, although to a lesser extent than the en¬ 
hancement seen for larger DM field voids. Another interest¬ 
ing aspect that is seen in the void abundances of the Galileon 
model is that the results of the full and linear variants are 
rather similar. This is very different from what is seen in the 
abundances of collapsed haloes, for which, due to the sup¬ 
pression effects of the screening mechanism, the full model 
has considerably fewer massive haloes than the linear variant 
(see e.g. Fig. 5 of Ref. IlhSll .I This illustrates that the effects 
of the screening mechanism are much weaker around under- 
dense regions, as expected. 

In the case of Nonlocal gravity, the number density of DM 
field voids is, within the etrorbars, the same in the full and 
QCDM variants. Here, recall that the largest voids found 
in the Nonlocal simulations are smaller than those in the 


Galileon simulations due to the smaller box size used. For 
instance, the largest DM field void found in the simulations of 
the Nonlocal model has R^ « 17Mpc//i. This, together with 
the fact that in the Galileon model the enhancement is most 
noticeable for R^ > 15Mpc/ft., suggests that the simulation 
box of the Nonlocal gravity model is not big enough to capture 
the impact of the fifth force on larger voids. Indeed, for large 
voids, there seems to be a trend for the full Nonlocal model 
to overpredict the number of voids with R^ > 15Mpc//i rel¬ 
ative to QCDM, although this is not significant due to the size 
of the errorbars. Nevertheless, for Ry ^ 15Mpc//i, the en¬ 
hancement in the full and linear variants of the Galileon model 
is already around ~ 10%, whereas in the Nonlocal model it 
is still consistent with zero. This shows that the effects of the 
modifications to gravity in the Nonlocal model are, in gen¬ 
eral, weaker than those in the Cubic Galileon, which is also 
expected. 

The results become noisier for voids found in the halo field 
due to the smaller number of tracers. For both the Galileon 
and Nonlocal gravity models, within the errorbars, the number 
density of voids is essentially the same in all model variants. 
However, as an exercise, if one ignores the size of the error- 
bars for a moment, then one notes that, at least qualitatively, 
the halo field voids show a similar behaviour to their DM field 
counterparts. In the case of the Galileon, for instance, the 
largest halo field voids, Ry > 40Mpc//i, are ^ 10 — 20% 
more abundant in the full and linear variants, compared to 
QCDM. This qualitative trend, backed up by the expectation 
based on physical intuition, suggests that with improved halo 
field void statistics one should recover, at least to a certain de¬ 
gree, the same physical behavior seen for the DM field voids. 

In the results that follow, we analyse our void catalogues by 
splitting them into two bins of radial size. We split the voids 
in the Galileon model according to 


DM field 

: bin 1 

DM field 

: bin 2 

Halo field 

: binl 

Halo field : 

: bin 2 


[5-12.5] Mpc/fi, 
[12.5 - 20] Mpc/h, 
[10 - 30] Mpc/fi, 
[30 - 50] Mpc/fi, 


Ry « 8.20Mpc//i, 
Ry « 14.2Mpc//i, 
Ry « 20.0Mpc/fi, 
Ry « 36.5Mpc/fi, 
(28) 


and in the Nonlocal model as 


DM field 

: bin 1 

DM field 

: bin 2 

Halo field 

: bin 1 

Halo field : 

: bin 2 


[4 — 10] Mpc/h, 
[10 — 16] Mpc/h, 
[10 - 20] Mpc/h, 
[20 - 30] Mpc/h, 


Ry « 7.30Mpc/fi, 
Ry « 11.4Mpc/fi, 
Ry « 14.5Mpc/fi, 
Ry 23.2Mpc/fi, 
(29) 


where Ry = Rv,i/Nhin is the mean void size in each bin, 
where Ry^i is the radius of the i-th void in the bin and IVbin 
is the number of voids in each bin. The exact value of Ry 
fluctuates only slightly (< 1%) in between the different model 
variants. We have found that this binning choice constitutes a 
good compromise between having enough voids in each bin, 
whilst making sure that the void properties do not vary too 
much within a bin. 
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FIG. 2. Cumulative void size function (number density of voids with radii above /?„) for the Cubic Galileon (left panels) and Nonlocal gravity 
(right panels) models. The upper panels show the number density of voids found in the DM (circles) and halo (squares) density fields for the 
full (blue), linear (green) and QCDM (red) variants of each model, as labelled. The lower panels show the difference relative to QCDM. The 
errorbars depict the variance across the five realizations of each variant. 


C. Void density profiles 

Figure shows the spherically averaged DM density field 
and halo density field void density prohles found in the sim¬ 
ulations of the three variants of the Galileon mode0 with the 
void sample split to the size bins according to Eq. ( [28] ). Figure 
l^is the same as Fig.[^ but for the Nonlocal model. The void 
density profiles are characterized by a density increase from 
R' = 0 towards i?' « 1; an overdense ridge at i?' ~ [1—1.5], 
which is associated with the filaments and walls that surround 
the void (the ridge is less pronounced for larger voids); and 
a steady decrease towards the cosmic mean, 5 = 0, at larger 
radii. In these hgures, the curves show the best-fitting profiles 
obtained using the five-parameter formula of Eq. ( [26] l, which 
fits the simulation results very well. Recall that the density 
profiles were computed using the DM density held (which is 
the mass density felt by photons), for both the voids identihed 
using the DM and halo distributions. This fact should be taken 
into account when comparing these prohles with others in the 
literature llTOl ITTI fT02l . 

The functional form of Eq. ( [26l l is inspired by the expres¬ 
sions proposed by earlier works II70l|II]| in the context of 


^ The average density profile of all the voids in each bin should be spherical 
to a good approximation, even though each individual void is not. 


ACDM. In particular, the formula proposed by Ref. ifTOll dif¬ 
fers from ours by hxing S2 = 1. This was used to ht to the 
density prohles of voids found from subsampled (i.e. diluted) 
DM tracer particle helds in ACDM (see Ref. fTOl for the de¬ 
tails). The authors of Ref. Gal further found that there are 
relations between the four free parameters of their formula, 
which can be used to effectively hx two of them. On the other 
hand, the formula proposed by Ref. Ga has S 2 = Si, and 
was used to ht the density prohles of voids constructed from 
mock and observed galaxy catalogues. In both of these works, 
the voids were found using watershed-based void hnders, as 
in this paper. Recently, Ref. nna explored the connection 
between the properties of voids found using watershed meth¬ 
ods and the predictions of theoretical models based on excur¬ 
sion set theory II 1001 . In Ref. 11021 . the authors also pointed 
out that the performance of the htting formulae proposed by 
Refs. llTOllTn may depend on some aspects of the analysis 
such as the tracer type, tracer density, dehnition of void cen¬ 
ter, etc. (see Ref. 11021 for the details). Compared to these 
other formulae, our void prohle of Eq. ( [26| may appear less 
appealing due to the fact that it has the extra free parameter 
S 2 - However, the flexibility that comes with S 2 is what allows 
our formula to be a very good ht to the simulation results, both 
for the voids found in the DM and halo helds, and for all the 
variants of the models of gravity we consider. Moreover, the 
differences relative to QCDM computed using the best-htting 
formulae also match very well the relative differences mea- 
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FIG. 3. Void density profiles, pmipm = 1 + 5, for the DM density field (circles, left panels) and halo density field (squares, right panels) 
voids found in the simulations of the three Galileon model variants (distinguished by the different colors, as labelled), plotted as function of the 
scaled radius R! = R/Rv The profiles are computed using the DM density field (which is the one felt by photons), for both the voids found 
in the DM and halo density fields. The solid (dashed) lines show the best-fitting density profiles, using the formula of Eq. j26| l, for the bin of 
smaller (larger) void sizes, as labelled. The bottom panels show the relative difference to QCDM. The errorbars depict the variance accross the 
five realizations of each variant. 
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FIG. 4. Same as Fig.|^but for the Nonlocal gravity model. 
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sured in the simulations (lower panels of Figs. |^and|g. We 
note that it is not the goal of this paper to determine if the 
voids found in our simulations are self-similar (i.e., indepen¬ 
dent of the void size) or universal (i.e., independent of tracer 
type and/or redshift). It is also not our goal to investigate the 
dependence of the void profiles on the value of the density at 
i? = 0 (in Ref. 01021 the void profiles are shown to depend 
also on this value, in addition to the dependence). From 
hereon, in our analysis, what is important to retain is that the 
void density profiles in the simulations are well described by 
Eq. (2^ which can therefore be used to calculate the force 
profiles and lensing signal. 

The impact of the fifth force is better seen when comparing 
the difference relative to QCDM. In the case of the Galileon 
model (Fig. [^, compared to QCDM, the voids in the full and 
linear variants are « 2 — 3% emptier in the inner regions, i.e. 
R' < 0.5, for both the DM and halo voids (although the re¬ 
sult is noisier for halo voids due to poorer statistics, specially 
for the smaller radius bin). Physically, this is because the en¬ 
hanced gravity favours the piling up of matter in the outer re¬ 
gions, leaving less matter inside the void. The fact that the pre¬ 
diction from the linear and full variants are so close illustrates, 
once again, that the effects of the screening mechanism are 
weak around voids. In Ref. BTll . similar results were found 
in the context of f{R) gravity, using a spherical underdensity 
based void finder 01031 . In particular, the authors of Ref. Il47l 
found that the voids in f{R) models can be up to « 5% emp¬ 
tier than in ACDM. In Fig. it is also worth noting that for 
the smaller radius bin of DM field voids, at R' ~ 0.5 — 1, 
the voids are more underdense in the linear variant than in the 
full model. We shall present an explanation for this in the next 
subsection, when we look at the force profiles in the Galileon 
model. 

The effects of the fifth force on the void profiles of the Non¬ 
local model (Fig.|^ are weaker than those seen in the Galileon 
case. In particular, for the DM field voids, the smaller void 
size bin in the full variant shows a decrement of only « 1%, 
relative to QCDM; the difference becomes consistent with 
zero for the larger size bin. In the case of the halo field voids, 
there is a systematic trend for the voids in the full Nonlocal 
model to be « 2 — 3% emptier than in QCDM for R' < 0.5, 
but the poorer statistics make it hard to draw any decisive con¬ 
clusions. Nevertheless, the result of Fig.j^shows that, overall, 
the fifth force effects on the void density profiles in the Non¬ 
local model are weaker than in the Galileon model, which is 
expected. 


D. Force profiles in the Galileon model 

Figure shows the force profiles of the voids in the vari¬ 
ants of the Galileon model. The circles, linked by the dotted 
lines, show the simulation results. These were obtained by 


spherically averaging the radial force field in the simulations, 
which was constructed by using the force information at the 
N-body particle positions. The solid curves show the analyti¬ 
cal result computed using the best-fitting void density profiles 
of Eq. ( |26] l (cf. Eig. [^. Eigure shows that, for the linear 
and QCDM variants, the analytical calculation is in very good 
agreement with the simulation results. However, the same is 
not true for the case of the full variant of the Galileon model. 
In this case, the analytical result differs from that of the sim¬ 
ulations for R' < 1.25, for DM held voids, and for R' < 1.0 
for halo held voids. More specihcally, for all cases shown, the 
analytical result of the full variant always predicts a stronger 
force (more negative) than the linear variant, inside the void. 
This result was already seen in Sec. |IIB[ when we analysed 
the behaviour of the regime (iii) discussed there (cf. Eig. [^. 
On the other hand, in the simulations, the force inside the 
smaller voids (left panels) of the full variant is weaker than in 
the linear case (i?' ~ 0.5 — 1). Eor larger voids (right panel), 
the full and linear variant simulations exhibit nearly the same 
force prohles. 

The reason why the forces in the simulations of the full 
model are weaker (less negative) than those computed analyt¬ 
ically using the spherically averaged density profiles can be 
linked to the effects of screening. The smooth void density 
profiles depicted in Eig. [^correspond only to an average den¬ 
sity field, which does not fully capture the detailed distribution 
of matter around the voids. A more realistic picture is that, in¬ 
side the voids and at their edges, there are higher density peaks 
associated with dark matter haloes and their respective infall 
regions. Close to these higher density regions, the fifth force 
in the Galileon model is suppressed by the screening mecha¬ 


nism (cf. regime (i) discussed in Sec. II Bi, which results in 
a weakening of the total forc^ Herein lies the explanation 
for the mismatch between the analytical result and the simula¬ 
tion force profiles. By averaging first the matter field, despite 
of the presence of higher density peaks, on average, one ends 
up with a smoother and lower density void profile. This pro¬ 
file, when used in the analytical calculation, gives a fifth force 
which is stronger in magnitude than the corresponding linear 
solution (cf regime (iii) in Sec. |II B| i. On the other hand, by 
averaging directly the forces at the particle positions, one is 
averaging a force field which is already affected by the sup¬ 
pression effects of the screening due to the existing higher 
density peaks. This is why the force profiles measured in the 
simulations are weaker (less negative) than the analytical re¬ 
sult, as seen in Eig.|^ In other words, since the force equation 
in the full model is nonlinear (cf Eqs. (pTj)), it makes a dif¬ 
ference whether one computes the force analytically from the 
averaged density field, or one computes the force by averaging 
directly the force field. In the case of the QCDM and linear 
variants, the force equation is linear, and as a result, the oper¬ 
ations of averaging the density and the force field commute, 
which is why there is almost perfect agreement between the 


® Even if one needs to fit the free parameters for different void sizes and for 
different density tracers. 


^ We note also that, in the simulations, the ad-hoc fix to keep the fifth force a 
real number is applied on a cell-by-cell basis on the adaptive mesh, which 
means that close to these density peaks the fix is not employed. 
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FIG. 5. Radial force profiles around the DM field (upper panels) and halo field (lower panels) voids in the variants of the Galileon model 
(distinguished by the different colors, as labelled). The circles with errorbars (which are in most cases smaller than the circles), linked by the 
dotted lines, correspond to the spherically averaged radial force field in the simulations. The solid lines correspond to the analytical prediction 
computed using the corresponding best-fitting void density profiles of Eq. shown in Fig.[^ The different panels show the result for the 
different void size bins, as labelled. What is actually plotted is the radial force scaled by the mean void size in each bin, T>,h /.R„ (cf. Eqs. 1 28 1 ). 
A negative sign for the force means that it points outwards. 


analytical and simulation results in these cases. 

Figure shows also that the suppression of the total force 
in the simulations of the full variant relative to the linear 
one is more pronounced in smaller voids. This is because 
smaller voids are denser, and therefore contain more higher 
density peaks per volume inside them and in their surround¬ 
ings, which enhances the suppresion effect of the screening. 
In particular, it is interesting to link this result with the dif¬ 
ferences between the linear and full model density profiles for 
the smaller size bin of the DM field haloes at R' ~ 0.5 — 1 
in Fig. As we noted in the previous section, on these ra¬ 
dial scales the voids in the linear model are slightly emptier. 
This can be explained by the fact that, in the simulations, the 
force in the linear model is stronger (more negative), which 
favours the evacuation of matter from inside the voicfl In 
principle, the same result should also be noticeable in the case 


It is worth noting that the relative differences in Fig. [^correspond to dif¬ 
ferent void populations, and as a result, some of the observed differences 
could arise from this. As a test, we have measured the density and force 
profiles in the full, linear and QCDM simulations, but at the spatial lo¬ 
cations of the voids in the QCDM model. This increases the chances of 
comparing voids that evolved from the same initial underdense regions. 
From this test we found only small quantitative changes with no impact on 
our conclusions. 


of the smaller size bin of the halo field voids, for which the 
force in the full model is also weaker than in the linear one 
at R' ~ 0.5 — 1 (lower left panel of Fig. [^. This is not visi¬ 
ble in Fig.|^(lower left panel), possibly because of the noisier 
measurements. 

For the Nonlocal gravity model, the force equations are lin¬ 
ear, as in the QCDM and linear variants of the Cubic Galileon 
model. As a result, there is good agreement between force 
profiles computed from simulations and the analytical results. 


IV. WEAK LENSING BY VOIDS IN MODIEIED GRAVITY 


In this section, we analyse the gravitational lensing signal 
from the voids in the Cubic Galileon and Nonlocal gravity 
models. We start by describing how to calculate the relevant 
lensing quantities in these theories of gravity, and then focus 
on the predictions using the voids found in the simulations of 
the two models. 
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A. Lensing shear calculation 

1. Lensing shear in GR 

The observable quantity in weak lensing studies II1041 1 1051 
is the reduced shear, g = 7t/(l — n), where 7t and k are 
called, respectively, the lensing tangential shear and conver¬ 
gence. The reduced shear is directly related to the elliptici- 
ties of the background galaxies whose light is distorted by in¬ 
tervening gravitational sources, which are voids in our case. 
In the weak-lensing regime, which is the regime for voids 
GMlIIOTl, 7t and k are both much smaller than unity, and 
consequently, one has g ve, 

The lensing convergence is obtained by integrating the 
Laplacian of the lensing potential <i)ien = (4) -f 'k) /2 as 

where I is the line of sight coordinate. Here, Eg = 
DgC^/(ATrGDdsDd) is called the critical surface mass density 
for lensing, where D^, Dg and Dds, are respectively, the angu¬ 
lar diameter distances between the observer and the void, the 
observer and the source galaxies, and the void and the source 
galaxies. We note, however, that given the way we choose to 
present our results below, the exact values of Ec are not im¬ 
portant (we comment further on this below). In GR, the two 
Newtonian potentials are the same T* = T* = $ien (in the 
absence of anisotropic stress). As a result, using the Poisson 
equation = AirCpm^, the convergence is given by 

J Sir',l')dl'= (31) 

where r' = r/Ry, I' = l/Ry, and r is a two-dimensional 
radial coordinate defined on the void plane (I = 0) with origin 
at the void center (i.e. R^ = r'^+P). From the above equation, 
one sees that in GR the lensing convergence is simply given 
by the projected density profile of the void, E(r) = i?„E(r'). 
We perform the integral of Eq. ([3T| l numerically, using the 
density contrast formula of Eq. \26). Note that the integral is 
dimensionless, and hence, k is dimensionless as well. 

The tangential shear is defined as 


')t = K - K,, 


(32) 


where 

2 r 2 r' 

/ yn{y)dy = Ry-^ / y'K{y')dy' (33) 
’’jo ’"Jo 

is the mean convergence inside radius r. The mean projected 
mass inside radius r is given by E(< r) = kEc. Eor clarity, 
we note that in the last equality of Eq. [^we have used that 
K{y) = Rytsiij'), where K{y') is the convergence given in 
terms of the radial coordinate scaled by Here, we follow 
Refs. Il49ll50ll and quote the lensing predictions in terms of the 
dijferential surface mass density AE(r), which is given by 

(34) 


By quoting the results in terms of AE(r) we avoid having to 
compute the values of E^. In Eqs. ( (3T] i and ( [3^ , we have also 
factored out the void radius, i?„, to make it explicit that the 
lensing quantities, depend linearly on it, e.g. AE oc Ry. 


2. Lensing in linear models of gravity with Ges(,z) 

As we have seen in Secs. |II A| and |IIB| in Nonlocal grav¬ 
ity and in the linear variant of the Cubic Galileon model, 
the gravitational law is the same as in GR, but with a scale- 
independent effective gravitational strength, GeS- Moreover, 
in these two models, the two Newtonian potentials are also 
the same and equal to the lensing potential, $ien = $ = T* 
(in the absence of anisotropic stress). Consequently, the dif¬ 
ferential surface mass density in these models is obtained in 
the same way as in GR, but taking the factor GeS/G properly 
into account. Explicitly, one has 

AE(r) = ^AE°^(r). (35) 

Gr 


3. Lensing in the full Cubic Galileon model 

As in GR, in the full variant of the Galileon model one also 
has that 4)ien = $ = 4*, and therefore, the lensing conver¬ 
gence is also given by integrating along the line of sight 


47rGE, 


-R, 


J (/,(')+ 2^(r',('))clT, 


(36) 


where /R and are given, respectively, by Eqs. (Ill 

and (121, and recall that in spherical coordinates, = 


+2$,rj /R. Given k, then the values of k and 7* are 
obtained as in Eqs. ( [3^ and ( [3^ , respectively. Here, we shall 
also quote the results for the full Galileon model in terms of 
the differential surface mass density, AE = ’Sc'jf However, 
one should bear in mind that the meaning of AE is different 
from the previous cases. In the full Galileon model, in ad¬ 
dition to the contribution from the projected mass, AE also 
depends on the projected distribution of the Galileon field 
(cf Eqs. ( [TT] ) and (|T^). 


B. Lensing by the voids in the simulations 

Figure shows the lensing signal associated with the halo 
field voids found in the simulations of the Galileon (left panel) 
and Nonlocal (right panel) gravity models. The curves were 
computed as described in the previous subsection using the 
best-fitting density profiles of Eq. ( [26| shown in the right pan¬ 
els of Eigs. |^and|^ Eor brevity, we show only the result 
for halo field voids. These are the ones that are more closely 
related to observations, where one first identifies voids using 
galaxy catalogues and then looks at the lensing signal at the 
void locations Il49ll50ll . Note also that the values of AE are 


AE(r) = E(< r) - E(r) = Ec7t. 
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FIG. 6. Lensing differential surface mass density, AE, for the best-fitting density profiles of the two void size bins (distinguished by the line 
styles, as labelled) for the halo density field voids found in the simulations of the different variants (distinguished by the colors, as labelled) of 
the Cubic Galileon (left panel) and Nonlocal (right panel) gravity models. The result is scaled by the inverse of the mean void radius in each 
size bin, . 


scaled by R~^, which means that the voids in the larger size 
bin (dashed curves) have lensing effects of larger magnitude 
(i.e. one needs to multiply the result by Ry). 

In the case of the Galileon model, the maximum of the lens¬ 
ing signal, which occurs at r' ~ 0.75 — 1, is approximately 
twice as strong in the full and linear variants, compared to 
QCDM. At r' < 1, the signal is slighlty stronger in the full 
than in the linear variant n This follows from the stronger 
fifth force (more negative) in the full model associated with 
the regime (iii) discussed in Sec. |IIB| However, in Fig.[^ we 
have seen that, due to the screening mechanism, the force pre¬ 
dictions of the full and linear variants measured in the simu¬ 
lations are actually closer than what is predicted by the ana¬ 
lytical calculation using the best-fitting density profiles. For 
this reason, it is reasonable to assume that the lensing signal 
of the full model is actually closer to the prediction of the 
linear variant, compared to what is observed in Fig. Never¬ 
theless, given that the differences between the full and linear 
variants are smaller than their differences relative to QCDM, 
it remains safe to conclude, as we have seen in previous sec¬ 
tions, that the effects of the fifth force in the Galileon model 
are quite pronounced in voids, where the screening is not very 
efficient. 

In the case of the Nonlocal model, the effects of the fifth 
force are considerably weaker than in the Galileon model. In 
particular, for the cases shown in Fig. the maximum ampli¬ 
tude in the value of |AE| (at r' ^ 0.75 — 1) is « 10% larger 
in the full Nonlocal model, compared to its QCDM variant. 
This illustrates that the effect of the modifications to gravity 
in the Nonlocal model are more challenging to detect using 
the lensing signal from voids. 

The modifications to gravity affect the lensing in voids in 
two main ways: (i) through the modifications to the average 


^ The "spikey" features seen in the blue curves follow from our ad hoc fix to 
the problem of complex fifth force solutions (cf. Sec. IffB). 


density profiles of voids; and (ii) directly through the modifi¬ 
cations to the lensing potential. Figure[7]measures the relative 
impact of these two effects in the Galileon (upper panel) and 
Nonlocal gravity (lower panel). In the figure, the red and blue 
curves are the same as in Fig. The black curves are com¬ 
puted using the lensing equations of the full variants of the 
Galileon and Nonlocal models, but using the density profile 
of the voids in the QCDM variants. As a result, comparing 
the red and black curves shows the effect of modifying the 
force law, whereas the difference between the black and blue 
lines shows the impact of the modified density profiles. In the 
case of the Galileon model. Fig. [7] shows that the dominant 
effect comes from the fifth force. This is seen by the large 
difference between the red and black curves. Figure]^ shows 
the result for the larger size bin of the halo field voids, which 
are slightly emptier in the full variant of the Galileon model, 
compared to QCDM (cf. bottom right panel in Fig. |^. This 
helps to further increase the amplitude of | AE|, but by a much 
smaller amount. On the other hand, for Nonlocal gravity, the 
direct effect of the fifth force on lensing is comparable to the 
effect of having slightly emptier voids (cf. Fig.|^. 

C. Connecting to observations 

References ll49l and ll50l have recently detected the lens¬ 
ing signal associated with voids in the galaxy distribution 
(see also Refs. Ill0614108ll for earlier forecast studies). In 
Ref 1491 . the authors stacked the voids of the catalogue of 
Ref ITOQI . which were found using a watershed algorithm in 
the three-dimensional main galaxy and luminous red galaxy 
(LRG) samples of the Sloan Digital Sky Survey-Data Release 
7 ms (SDSS-DR7). On the other hand, in Ref iH, the au¬ 
thors used also the SDSS LRG catalogues, but the voids were 
found using a method that is specifically designed for lens¬ 
ing. In this method, emptier regions are found in projected 
two-dimensional slices of the survey volume, which seems to 
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FIG. 7. Relative impact of the fifth force and modified density pro¬ 
files on the lensing signal of halo field voids in the Cubic Galileon 
(upper panel) and Nonlocal (lower panel) gravity models. For both 
models, the red and blue lines have the same meaning as those in 
Fig. § The black curves are obtained by calculating the lensing sig¬ 
nal with the full Galileon and full Nonlocal model force equations, 
but using the best-fitting void density profile of the QCDM voids. 
The comparison between the red and black measures the effects of 
the fifth force alone; whereas comparing the black and blue curves 
shows the impact of the modified void density profiles. 


increase the significance of the lensing detection. 

As we discuss below, a robust comparison between these 
observations and the results of Fig. [^requires more detailed 
modelling of the theoretical predictions. Nevertheless, one 
can still compare some of our results to try to get a feeling 
about what these measurements imply for modified gravity. 
For instance, in Fig. 5 of Ref. mni, it is shown that for voids 
with size G [15,30] Mpc/h, the values of the differen¬ 
tial surface mass density at its minimum are, approximately, 
within AS S [—0.4, —0.7] lO^^MQh/Mpc^ (this estimate is 
based on the size of the errorbars there). From Fig. [^ for the 
case of the smaller size bin of the full Galileon model we have 
min (AS) « -0.065i?„ = -1.3 x 10^^MQh/Mpc^._ For 
the full Nonlocal model, we have min (AS) « —0.032i?„ = 
-0.46 X lO^^Afgh/Mpc^ and min (AS) « -0.039i?^ = 
—0.9 X IO^^Mq/i/Mpc^, for the smaller and larger size bins, 


respectively. Hence, for both the Galileon and Nonlocal mod¬ 
els of gravity, we get the same typical order of magnitude as 
in the observations. One notes that in the case of the Galileon 
model, the size of the effect is larger than the results presented 
in Ref. ifSOll . This suggests that, indeed, lensing by voids may 
have the potential to constrain models like the Galileon. 

Before summarizing our results in the next section, we find 
it instructive to briefly comment on a number of aspects that 
should to be taken into account before properly confronting 
these (and other) models to lensing observations. These as¬ 
pects include; 

1. Impact of Sc In Fig.j^ we quote our results in terms of AS, 
but in reality, what one measures directly from galaxy el- 
lipticities is the shear, p « 7^ = AS/Sc. The calcula¬ 
tion of Sc depends on the cosmological background, which 
can be different between the Galileon, Nonlocal, and the 
standard ACDM models. Consequently, if in observational 
studies, one measures 74, but quotes the results in terms of 
AS by assuming a background cosmology to compute Sc, 
then this may introduce some bias that should be carefully 
addressed. Furthermore, Sc depends also on the redshift 
distribution of the source galaxy population, although this 
can always be set accordingly using the properties of the 
observed galaxies. 

2. Void redshift distribution The lensing signal in Fig. was 
obtained analytically using the density profiles of the voids 
in the simulations at z = 0. In the observations, however, 
the lensing signal is detected by stacking voids that span a 
given redshift distribution z > 0. In the particular case of 
the Galileon and Nonlocal gravity models, the fifth force 
is weaker at earlier times (see e.g. Fig. 3 of Ref. Il68ll and 
Fig. 2 fo Ref. l56l ). which reduces the amplitude of the 
signal depicted in Fig. [^ 

3. Void stacking The lensing signal associated with individ¬ 
ual voids is too weak to be detected in current observa¬ 
tions, which is why Refs. iilEol used stacked voids in 
their analyses. When interpreting such results in modified 
gravity, for a given stack, voids at different redshifts have 
different weights in the observed lensing signal because of 
the redshift dependence of the fifth force, Sc and also of 
the screening efficiency. Such effects should be taken into 
account if one, for instance, tries to use the lensing obser¬ 
vations to reconstruct a mean density profile for the stack. 
Here, an interesting analysis would involve stacks of voids 
binned not only by size, but also by redshift. 

4. Systematic biases The lensing calculations performed here 
assume that the density distribution in voids is perfectly 
smooth. In reality, however, voids contain substructure and 
its amount is expected to be different in models with differ¬ 
ent growth rates of structure. Given that the lensing signal 
from voids is relatively weak (compared to that induced by 
DM haloes) it may be interesting to investigate the extent to 
which void substructure can impact on the overall lensing 
signal. This can be studied by looking at the lensing sig¬ 
nal using ray-tracing methods in the simulations, without 
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modelling their profiles as a smooth distribution. Our lens- 
ing calculations also assume that the void is the only source 
of lensing. A ray-tracing analysis would also help to bet¬ 
ter quantify the contamination of the lensing signal coming 
from intervening matter along the line of sight. 


5. Screening effects Related to the above point, a ray-tracing 
analysis is also able to capture more accurately the effects 

we 


of the nonlinear screening mechanism. In Sec. HID 


saw that, in the full variant of the Galileon model, it makes 
a difference whether one computes the force profiles analyt¬ 
ically from the spherically averaged density profiles, or by 
spherically averaging the force field directly. Moreover, the 
efficiency of the Vainshtein mechanism depends also on the 
geometry of the mass distribution as investigated recently 
in Refs. Il23ll24l . This means that calculations based on the 
mean spherical profile of a stack of voids may not fully cap¬ 
ture the fifth force effects from each individual nonspheri- 
cal void. These issues can be circunvented by directly in¬ 
tegrating along the line of sight for each void using 
ray-tracing and stacking the resulting signal. In this way, 
one probes directly the lensing potential distribution with¬ 
out introducing any bias that arises when one averages first 
the density field. 


6. Halo/galaxy abundance and bias We have found voids us¬ 
ing both DM particles and DM haloes as tracers, with the 
latter case being that which is more relevant when compar¬ 
ing to observations. Due to halo bias, haloes cluster differ¬ 
ently depending on their mass, and hence, the resulting void 
catalogues depend on the minimum halo mass cut used to 
identify them. This is turn has an impact on the abundances 
and profiles of voids BtI 111 II 11121 . Since different types 
of galaxies populate haloes differently, a robust comparison 
with observations should ensure that the number density 
and bias of the tracers used in simulations matches those 
of the tracers used in observations. A first step towards this 
goal could be to construct mock galaxy catalogues using 
halo occupation distribution modelling II113I . as Ref. Una 
has done for ACDM. Such an analysis can also tell whether 
voids identified using certain types of density tracers are 
better suited for tests of modified gravity. 


7. Combining different void finders The way voids are found 
in simulations and/or in real galaxy catalogues can also 
affect the resulting lensing signal. For instance, as we 
mentioned above, the authors of Ref. Eol optimize their 
analysis for lensing by finding the voids in projected two- 
dimensional slices of a spectroscopic galaxy survey. This 
may partly explain why the significance of their detection 
is higher than that found in Ref. 1491 . in which the voids 
are found in three dimensions. It would therefore be of 
great interest to find voids in the way of Ref. Eol in the 
N-body simulations as well. It is also well known that 
different void finding techniques yield different void pro¬ 
files iTon . and hence, different lensing predictions. For in¬ 
stance, voids found with spherical underdensity (SU) meth¬ 
ods have sharper transitions from the inside of the void to 
the surrounding ridge. This boosts the lensing effect, as can 


be checked in the a panel of Fig.j^in the Appendix, where 
a is the parameter of the formula of Eq. ( [26] l that controls 
the slope of this transition. Moreover, the recent work of 
Ref. ll 14l has shown that it may be more natural to char¬ 
acterize the void profiles with respect to their boundaries 
(which is where most of the mass is), instead of with re¬ 
spect to the void center (which is devoid of tracers). This 
also results in steeper density profiles close to the void edge 
(see Fig. 7 of Ref. 01141 1. The differences in void profiles 
obtained with different void finding methods is generally 
portrayed as a source of uncertainty in void related work, 
but we note that some advantages may arise from it. For in¬ 
stance, since the fifth force acts to make voids emptier and 
the ridges denser, then methods like SU or that of Ref. 01140 
may be particularly suitable for modified gravity studies, 
as they may amplify the size of the fifth force effects (see 
e.g. Ref. Ell, where the authors use SU methods to study 
voids in f{R) gravity). Hence, we believe that the com¬ 
bination of the results from different void finding methods 
(provided they are consistently used in simulated and real 
data) is something to be explored with more detail when 
designing observational tests. These investigations are the 
subject of ongoing work ma. 


V. SUMMARY & CONCLUSIONS 


We have studied the lensing signal associated with voids in 
Cubic Galileon and Nonlocal gravity cosmologies, which are 
examples of models that modify the gravitational lensing po¬ 
tential. The gravitational law in the Nonlocal model can be 
parametrized by an enhanced effective gravitational strength 
(Geff « 1.06G, at z = 0), which is independent of the length 
scale. In the Galileon model, the modifications to gravity are 
scale-dependent and in Sec IIB we discussed three relevant 
regimes: (i) a regime which occurs close to massive bodies, in 
which the fifth force is suppressed via the Vainshtein screen¬ 
ing mechanism; (ii) a regime which occurs in regions of small 
density constrast, |<5| <C 1, where the equations become linear 
and the total force can be parametrized by an effective grav¬ 
itational strength, (GeS ~ 2G, at z = 0); and finally, (iii) a 
regime which occurs in regions where the density contrast be¬ 
comes sufficiently negative, where the amplitude of the fifth 
force is the largest (more negative). 

The fifth force in the Galileon and Nonlocal gravity mod¬ 
els has an impact on the lensing signal in and around voids 
through two main effects. First, the fifth force changes the 
density profiles of the voids, and second, it also modifies the 
lensing potential directly. This means that even for fixed mass 
distribution, the lensing signal in these theories of gravity is 
still modified w.r.t. GR. This is different from other popular 
models like f{R) and/or DGP gravity, which practically do 
not directly modify the lensing potential. Hence, models that 
directly modify lensing are more amenable to being tested by 
lensing observations than those that do not. 

We have used results from N-body simulations to study the 
effect of the fifth force in these two theories. We analysed 
the abundances and profiles of the density, force and lensing 
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shear of the voids found in the DM and halo density fields of 
the simulations using a watershed algorithm. When assessing 
the impact of the modifications to gravity in these two mod¬ 
els, we always compared their predictions to models called 
QCDM, which have the same background expansion as the 
respective Galileon and Nonlocal models, but have GR as the 
theory of gravity. In the case of the Galileon model, we have 
also analysed the results of a model variant with a linearized 
scale-independent force law (cf. Tablej^. Our main results can 
be summarized as follows; 

• In the Galileon model, the fifth force boosts the abun¬ 
dance of the larger radius DM field voids {R^ > 15Mpc//i) 
by « 10% — 30% (cf. Fig. [^. This is because the enhanced 
gravity causes voids to expand faster and also favours the 
merging of smaller voids into larger ones. For the voids found 
in the halo density field, the same qualitative trend is also seen 
but is less pronounced due to poorer statistics. In the case 
of the Nonlocal model, the modifications to gravity are not 
strong enough to leave a clear signal on the abundances of the 
voids found in our simulations (cf. Fig.[^. 

• Our five-parameter formula of Eq. ( [26l l fits very well 
the DM and halo field void density profiles in the simulations 
for all the model variants and for a wide range of void sizes 
(cf. Figs. 1^ and |^. Our formula contains an extra parame¬ 
ter compared to others used recently in the literature GQlliI], 
which gives it the flexibility required to provide good fits. We 
used our best-fitting formula to compute analytically the lens- 
ing signal associated with the voids in the simulations. In¬ 
vestigations about the self-similar or universal nature of the 
profiles in our simulations are left for future work. 

• The fifth force in the Galileon model makes the voids 
slightly emptier (« 2% — 3%) in their inner parts, R' < 
0.5 (cf. Fig. s- This is because the enhanced gravitational 
strength favours the evacuation of matter from the inside of the 
void into the surrounding filament and wall structures. This 
result is seen for both DM and halo field voids, although the 
signal is more significant for the smaller DM field voids. In 
the case of the Nonlocal model, the gravitational strength is 
also enhanced, and so one expects the same qualitative be¬ 
havior. Quantitatively though, the weaker fifth force in this 
model, together with the size of the errorbars allowed by our 
simulations, makes it more difficult to see the effects of the 
modifications to gravity (cf. Fig.|^. 

• Inside the voids of the full Galileon model, the force 
measured directly from the simulations is weaker (less nega¬ 
tive) than the force computed analytically from the best-fitting 
void density profile (cf. Fig.|^. We have attributed this to the 
screening by high density peaks that exist inside the voids and 
in their surroundings, whose effect gets diluted if one averages 
the density field first to compute the force analytically. On the 
other hand, for the linear and QCDM variants, the analytical 
result is in very good agreement with the force measured di¬ 
rectly from the simulations. 

• The effects of the fifth force in the Galileon model 
can make the lensing signal in voids approximately twice as 
strong as in GR (cf Fig.|^. This large difference comes pre¬ 
dominantly from the modifications of the lensing potential per 
se, with the different void density profiles being of secondary 


importance (cf. Fig. |^. In the case of the Nonlocal gravity 
model, the fifth force also enhances the expected suppression 
of the lensing signal, but only by « 10% (cf. Fig.[^. In this 
model, the modifications to the density profiles and direct ef¬ 
fects of the fifth force on lensing contribute equally to the dif¬ 
ference relative to GR (cf. Fig.|7]). 

• For all our Galileon model results, the predictions from 
the full and linear variants are of comparable size. This is dif¬ 
ferent from the case of predictions associated with dark matter 
haloes (like their abundances or concentration), for which the 
effects of the full variant are typically much smaller than those 
of the linear variant because of the screening 1371 . This illus¬ 
trates that the suppression effects of the screening mechanism 
are not very strong around voids, which is why the latter can 
be regarded as potentially powerful probes of gravity on cos¬ 
mological scales. 


Overall, the results in this paper show that observations of 
the lensing signal associated with voids can prove very valu¬ 
able in constraining the gravitational law on large scales. In 
the future, we plan to use some of the results presented here to 
help to develop more robust observational tests, by following 
the steps outlined in Sec. IV C We believe that such investi¬ 
gations would be timely, specially when interpreted in light 
of future observational missions such as DESI Ifira . LSST 
ifTTTi and Euclid ifTTSl . 
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Appendix A: Parameter impact in the void fitting formula 


Figure shows the effect that each of the five parameters 
{Sv,a, (3, si, S 2 ) that enter Eq. (261 have on the void density 
and void lensing differential mass density profiles. The cal¬ 
culation of the lensing signal was performed with GR as the 
theory of gravity. 
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FIG. 8. Impact of each of the parameters of the void density contrast formula of Eq. \2()) (upper panels). The bottom panels show the respective 
lensing signal AS. The curves are colored by the values of the parameter that is varying in each panel (from left to right, these are, 5^, a, P, 
Si and S 2 , respectively). The parameter values are indicated by the color bar at the right of each panel. When one parameter varies, the others 
are held fixed at their base values which are a, P, Si, S2) = (—0.6,3, 7, 0.9,1.1). The calculation of the lensing signal was performed 
assuming GR. 
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